%% Series 6 - Exercise 1 for the course of DOE
% Sensitivity of a Placket-Burman design
%
% Table of Contents
% Statement
% Establishing the base of the exercise
% Linear model without interactions
%     Analysis
% Linear model with interactions
%     Observations
% Test with the full fold-over
%     Observations
% The function simulating the measurement
%% Statement
%
% To analyse the sensitivity of a Plackett-Burman design, the measurement of a linear system with first-order interactions is simulated and the impact of interactions and uncertainties on the main effects is observed.
% Start by coding a  Matlab function, Y=measurement(E,[,sigma,alpha]), simulating a series of experiments whose true model is
%         y=a_o+\sum_{i=1}^{5} a_i x_i + \sum_{i \leq j}^{5} a_{ik} x_j x_i
% and whose experimental uncertainty is modeled by the distribution Normal distribution N(\mu=0,\sigma^2). The parameters of this function must be the matrix of essays E and the standard error \sigma. The coefficients \alpha of the true model can be given as parameters of the function or defined by default within the function (use the logical variable nargin).
% After that, in the main part of the script :
% prepare a Plackett-Burman design,
% simulate the essays with the routine measurement(),
% analyze the results to determine the effets with the routine fitlm().
% c) Perform this analysis for a true model without interaction with important (10%), medium (5%) and small effects (1%) varying the uncertainty of the measurement between 1% and 10% of the mean results. Observe the evolution of the estimated coefficients, of the confidence intervals and of the p-values of the ANOVA.
% d) Perform the analysis for a true model whose levels of interaction vary and with the measurement uncertainty varying between 1% and 10%. Observe the evolution of the estimated coefficients and confront these results with the matrix of alias between the main effects and the interactions.
% e) Finally perform a full foldover to get rid of the influence of the interactions.
%
%% Establishing the base of the exercise
% The function measurement() is coded at the end of this file
% Defining a true model without interaction purely linear whose coefficients are

coef=[100;
    -5; -10; 12; -1; 2; ...
    0; 0; 0; 0; ...
    0; 0; 0; ...
    0; 0; ...
    0 ];

coefNames={'a_0'   'a_1'   'a_2'   'a_3'   'a_4'   'a_5'};
disp(array2table(coef(1:6)',"VariableNames",coefNames))

% Building an Hadamard design H_8

X=hadamard(8);
varNames={'I'   'x1'   'x2'   'x3'   'x4'   'x5'   'x6'   'x7'};
disp(array2table(X,"VariableNames",varNames))

% Selecting a matrix of experiment corresponding to the last 5 columns of the Hadamard matrix

E=X(:,2:6);
disp(array2table(E,"VariableNames",varNames(2:6)))

% This experimental situation and the set of experiments is now defined. Let's vary the experimental parameters to observe the sensitivity of the Plackett-Burman design.
%%
%% Linear model without interactions
% Varying the standard deviation of the measurements between 1% and 9%

for i=1:2:10
    sigma=i; % Choose the standard deviation
    Y=mesure(E,sigma,coef); % Run the experiments
    
    
    mdl_lin=fitlm(E,Y,'linear'); % Analyze the results (linear model)
    
    Effects=mdl_lin.Coefficients.Estimate;
    
    figure
    subplot(211)
    plotEffects(mdl_lin)
    title(['Effects, s= ' num2str(sigma)])
    
    tbl=anova(mdl_lin,'component',2);
    disp(' ')
    disp(['ANOVA (SStype=2) for s=' num2str(sigma)])
    disp(' ')
    disp(tbl)
    
    CI=array2table(coefCI(mdl_lin),...
        'RowNames',coefNames(1:6),...
        'VariableNames',{'min' ,'max'});
    
    disp(' ')
    disp(['Confidence intervals of the coefficients for s=' num2str(sigma)])
    disp(' ')
    
    disp(CI)
    disp(' ')
    
    subplot(212)
    bar([coef(2:6) Effects(2:end) ])
    legend('True', 'Estimated')
    title(['Half-effets, s= ' num2str(sigma)])
    
end

%% Analysis
% You may observe that the accuracy of the estimation is decreasing with the increase of the standard deviation. When the standard error is 1% of the mean response, all  coefficients are well determined excepted a_4 whose magnitude is equivalent at the magnitude of the standard error.
%  The estimation of the coefficients stays stable even when the variance has increased enormously.
% It is worth to run the code several times to observe the changes due to the experimental error .
%
%% Linear model with interactions
% A true model with interactions
% Variance varying between 1% and 9%
%               The coefficients of the true model are:

coef=[100;
    -5; -10; 12; -13; 2; ...
    10; 5; 0; 0; ...
    0; 0; 0; ...
    0; 0; ...
    0 ];

coefNames={'a_0'   'a_1'   'a_2'   'a_3'   'a_4'   'a_5' 'a_12'   'a_13'};
disp(array2table(coef(1:8)',"VariableNames",coefNames))

% Building an Hadamard design H_8

X=hadamard(8);
E=X(:,2:6); % Select the matrix of experiment

for i=1:2:10
    sigma=i; % select the standard deviation
    Y=mesure(E,sigma,coef); % simulating the experiments
    
    mdl_lin=fitlm(E,Y,'linear'); % analysing the results
    
    Effects=mdl_lin.Coefficients.Estimate;
    figure
    subplot(211)
    plotEffects(mdl_lin)
    title(['Effects, s= ' num2str(sigma)])
    tbl=anova(mdl_lin,'component',2);
    disp(' ')
    disp(['ANOVA (SStype=2) for s=' num2str(sigma)])
    disp(' ')
    disp(tbl)
    
       CI=array2table(coefCI(mdl_lin),...
        'RowNames',coefNames(1:6),...
        'VariableNames',{'min' ,'max'});
    disp(' ')
    disp(['Intervals of confidence of the coefficients for s=' num2str(sigma)])
    disp(' ')
    disp(CI)
    disp(' ')
    subplot(212)
    bar([coef(2:6) Effects(2:end) ])
    legend('True', 'Estimated')
    title(['Half-effects, s= ' num2str(sigma)])
    
end

%% Observations
% You may observe that independently from the variance, the results are systematically biased by the interactions that are aliased with some of the coefficients .

% model matrix for the two non-zero interactions
X2=x2fx(E,[1 1 0 0 0;1 0 1 0 0]);
A=array2table(X'*X2/8,...
    "RowNames",{'lo' 'l1' 'l2' 'l3' 'l4' 'l5' 'l6' 'l7'}, ...
    "VariableNames",{'a12' 'a13'});
disp(A)

%  The main effects a_1, a_4 and a_5 which are not aliased behave as in the pure linear case.
% The contrasts l_6 and l_7 corresponding to dummy variables,  by a smart re-ordering of the factors could be used as estimators of  the two interactions, if they are expected. But this is away of the objective of this exercise.
%
%% Test with the full fold-over
% The true model has interactions
% The variance change succesively from 1% to 9%
% The design is a full foldover
%                The coefficients of the true model are:

coef=[100;
    -5; -10; 12; -13; 2; ...
    10; 5; 0; 0; ...
    0; 0; 0; ...
    0; 0; ...
    0 ];

coefNames={'a_0'   'a_1'   'a_2'   'a_3'   'a_4'   'a_5' 'a_12'   'a_13'};
disp(array2table(coef(1:8)',"VariableNames",coefNames))
X=hadamard(8);
E=X(:,2:6);
E= [E;-E]; % full foldover
for i=1:2:10
    sigma=i; % determine the standard deviation
    Y=mesure(E,sigma,coef);
    mdl_lin=fitlm(E,Y,'linear');
    Effects=mdl_lin.Coefficients.Estimate;
    figure
    subplot(211)
    plotEffects(mdl_lin)
    title(['Effets, s= ' num2str(sigma)])
    
    tbl=anova(mdl_lin,'component',2);
    disp(' ')
    disp(['ANOVA (SStype=2) for s=' num2str(sigma)])
    disp(' ')
    disp(tbl)
    
     CI=array2table(coefCI(mdl_lin),...
    'RowNames',coefNames(1:6),...
    'VariableNames',{'min' ,'max'});
    disp(' ')
    disp(['Intervals of confidence of the coefficients for s=' num2str(sigma)])
    disp(' ')
    disp(CI)
    disp(' ')
    subplot(212)
    bar([coef(2:6) Effects(2:end) ])
    legend('True', 'Estimated')
    title(['Demi-effets, s= ' num2str(sigma)])
    
end

%% Observations
% With to the full foldover design the problem due to the interactions has disappeared. The Plackett-Burman design is robust against a sample variance increase, but not against the possible presence of interactions. The latter must be mitigate with the full fold-over technique.
%
%% The function simulating the measurement

function Y=mesure(E,sigma,coef)

% E matrix of experiments with 5 factors
% Y experimental results
if nargin==1 % default variance
    sigma=2;
end

if (nargin==2) || (nargin==1) % default coefficents
    coef=[100;
        8; -10; 12; -1; 2; ...
        0; 0; 0; 0; ...
        0; 0; 0; ...
        0; 0; ...
        0 ];
end
spec=[ 0 0 0 0 0
    1 0 0 0 0
    0 1 0 0 0
    0 0 1 0 0
    0 0 0 1 0
    0 0 0 0 1
    1 1 0 0 0
    1 0 1 0 0
    1 0 0 1 0
    1 0 0 0 1
    0 1 1 0 0
    0 1 0 1 0
    0 1 0 0 1
    0 0 1 1 0
    0 0 1 0 1
    0 0 0 1 1];

% Model matrix
X=x2fx(E,spec);


% Generating the error
er=randn(size(E(:,1)))*sigma;

% Generating the results
Y= X*coef + er;

end